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Abstract. 

We present preliminary results from self-consistent, high resolution direct A^-body 
simulations of massive black hole binaries in mergers of galactic nuclei. The dynamics of the 
black hole binary includes the full Post-Newtonian corrections (up to 2.5PN) to its equations of 
motion. We show that massive black holes starting at separations of 100 pc can evolve down to 
gravitational-wave-induced coalescence in less than a Hubble time. The binaries, in our models, 
often form with very high eccentricity and, as a result, reach separations of 50 Schwarzschild 
radius with eccentricities which are clearly distinct from zero — even though gravitational wave 
emission damps the eccentricity during the inspiral. These deviations from exact circular orbits, 
at such small separations, may have important consequences for LISA data analysis. 



1. Introduction 

Massive black holes (henceforth MBHs) are ubiquituous in the centers of galaxies. During the last 
decade, a great deal of effort was made observationally to characterize the physical properties of 
galactic nuclei and their black holes [1]. A number of phenomenological relations were unveiled: 
the masses of the central MBH appear to be correlated with some physical properties of the host 
nucleus, namely the luminosity [2] and the velocity dispersion [21 H]. Massive galaxies, however, 
typically do not evolve in isolation and undergo several merger events during their lifetime. If 
each of the merging galaxies harbors a MBH at the center, then a MBH binary is expected to 
form in the core of the resulting remnant. 

The paradigm for MBH binary evolution consists essentially of three distinct phases [5j : (i) 
Right after the merger the two MBHs sink on a dynamical friction time scale, towards the 
center of the remnant where they form a bound pair; (ii) Once such bound pair is formed, 
its semimajor axis continues to shrink by dynamical friction until the binary becomes hard; 
afterwards, the binary continues to harden by the ejection of field stars; (iii) When the binary's 
separation becomes small enough, gravitational wave (henceforth GW) emission may drive the 
final stage of relativistic inspiral and coalescence. The transition from (i) to (ii) is understood 



to be unproblematic, as it typically occurs on a cosmologically short time scale (< 10^^^ yr.) 
as long as the mass ratio is q = M1/M2 > 0.1 [6]; the subsequent transition from Phase (ii) to 
Phase (iii), however, could constitute a bottleneck for the evolution towards relativistic inspiral 
and coalescence. The binary stalling problem is especially severe in the case of both purely 
stellar (ie. with no gas), spherically symmetric galaxy models. In this case, the pool of stars 
available for close interaction with the MBH binary is very rapidly exhausted (roughly on the 
order of a local dynamical time scale) , and leaves the binary dynamically isolated at separations 
of the order of a parsec, corresponding roughly to the "hard binary separation" = G fired/ ^'^'^ , 
where fired is the binary's reduced mass and a is the velocity dispersion of the surrounding stellar 
field. This is the so-called Final Parsec Problem (hereafter FPP). If the time scale on which 
stars can diffuse in angular momentum space to enter the binary's loss cone is longer than a 
Hubble time, then the MBH binary evolution will stall and it will not become a GW source 
detectable by LISA. There are many proposed solutions to the FPP: dissipation of the binary's 
orbital energy in a gaseous circumbinary disk [318], massive perturbers [9], the infall of a third 
(Intermediate-) MBH, a triaxial potential may keep the loss cone full \10 \ [T2]. etc. 

Here we present preliminary results of A'-body simulations of idealized models of merging 
spherical galactic nuclei, starting from an initial separation of tens of parsec. We follow the 
subsequent formation and evolution of the MBH binary until it reaches the final stages of 
relativistic inspiral at separations of a few Schwarschild radii (~ few xlO"'' pc, if the MBHs 
have masses ~ IO^Mq). The nucleus that results from the merger has a pronounced triaxial 
structure; triaxiality facilitates the MBH binary's transition to the GW regime, presumably by 
inducing a collisionless regime which is capable of repopulating the loss cone due to centrophilic 
orbits. In sec. 2, we describe the galaxy merger set-up and the formation of the bound MBH 
pair. In sec. 3, we show that the MBH binary reaches the GW-dominated inspiral phase while 
the merger remnant still retains some degree of triaxiality. In sec 4, we show that the MBH 
binary may reach the LISA band with an eccentricity different from zero. We conclude with 
Sec. 5. 

2. Simulation of the merger of two spherical galactic nuclei 

A number of studies were dedicated to the FPP for MBH binaries — placed in the center of 
a spherical model of a galaxy merger remnant — using the Fokker-Planck formalism, A'-body 
simulations, or both [131 [HI [151 [HI [T71 US] . In those models, the loss cone is mostly empty and 
repopulation of loss cone orbits is driven by two-body relaxation in the diffusive regime [19\ I14j. 
The characteristic time scale for this process is the relaxation time, which is proportional to the 
total number N of stars in the nucleus. For the number of stars present in a typical galactic 
nucleus, this generally exceeds the Hubble time by a large margin — the exception being the 
most compact and dense nuclei. From cosmological studies, we expect that the mass function 
of MBH binaries peaks in the 10^ - W^Mq range and at high redshift z > 4 [20]; in such cases, 
the diffusive refilling of the loss cone may be marginally efficient, even in the spherical case, to 
drive the MBHs to coalescence in less than a Hubble time [18j. 

Nevertheless, the spherical approximation is not only a worst case scenario from the point 
of view of the loss cone dynamics, but it is also highly unrealistic in the sense that no merger 
remnant will ever exhibit such a high degree of symmetry. It is known that some significant 
(transient) triaxiality results almost inevitably from a galaxy merger [21] — as well as some net 
rotation, details depending on stellar and gas content, orbital parameters of the merger, etc. We 
have recently modelled the galactic nucleus remnant with rotating King models; in some cases, 
the rotation was high enough that the model was unstable to the formation of a rotating, triaxial 
bar- like feature. In such setting, the evolution of the MBH binary did not stall at a ~ a/i, but 
did inspiral down to coalescence on timescales well below a Hubble time [11^ [T2] . 

Notwithstanding the improvement of the latter models as compared to spherical ones, we 



have decided to take the next step further and set up a number of galaxy mergers p2]. We 
let two galactic nuclei to merge on nearly parabolic orbits, since this seems to be a typical 
configuration for mergers as is observed in state-of-the-art cosmological simulations [23]. Each 
galactic nucleus is represented by a spherically symmetric model with a power law density cusp 
of logarithmic slope 7 at the center, and a break radius which we set to unity in model units 
[24^ I25j: incidentally, this also permits us to study the mass deficits resulting from the ejection 
of stars by the binary. A massive point particle represents the MBH which, at the beginning, is 
placed in the center of each nucleus with zero velocity with respect to each nucleus. The total 
mass in stars of each galactic nucleus is unity; we adopt units where G = 1. In our sample, we 
have currently some sixty simulations in which we vary essentially four parameters: the total 
mass of the MBH binary M12 = Mi + M2 in units of the total stellar mass in each merger 
progenitor, its mass ratio q = M1/M2 < 1, the central slope 7 of the stellar density of each 
galaxy and the total number = A^i -|- A''2 of stars in the two galaxies. We have generated a 
sample of simulations for the following parameter range: 5 x 10~^ < M12 < 6 x 10~^, 0.2 < q < 1, 
0.5 < 7 < 1.2 and QAK <N< 256K. 
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Figure 1. Histogram representing the 
eccentricity distribution with which the MBH 
binaries settle when they form a bound pair. 
This value is an average value taken over a few 
A-body time units, after the binary becomes 
"hard". 



Figure 2. Fractional value of the per- 
turbing accelerations, normalized by the 
Newtonian acceleration between the MBHs. 
From: (a) the stars (green), (b) the to- 
tal lPN-h2PN-F2.5PN terms (blue), (c) the 
2.5PN term only (magenta). The perturbing 
forces are always orders of magnitude weaker 
than the dominant Keplerian force. 



After the two nuclei have merged, the MBHs sink to the center of the remnant where they 
quickly form a bound pair. In Fig. 1, we show the statistics for the average eccentricity 
with which the MBHs form a bound pair in our sample. There is a broad range of possible 
eccentricities, but with a clear bias towards the formation of binaries with very high eccentricity 
e > 0.8. This tendency for the formation of highly eccentric binaries results naturally from the 
(cosmologically motivated) initial conditions with which we set up the merger: the nuclei fall into 
each other on near-parabolic orbits; therefore the MBHs, being much heavier than the stars, also 
approach each other on near-parabolic orbits despite the perturbations from the surrounding 
stellar field. In our previous work, we have modelled the galactic nucleus of a merger remnant 
with rotating King models; in such models, the MBH binaries also tend to become bound with 
very high eccentricity I12j . In the latter work, it has been shown that the hardening rate 
due to the superelastic scattering of field stars is essentially independent of the eccentricity. It 
is well known that the emission of GWs has a very steep dependence on the pericenter distance 



and thus a highly eccentric binary inspirals on a much shorter time scale [26l 



3. Orbital evolution of the massive black hole binary 

We have implemented the relativistic effects to the MBH binary only, by using the Post- 
Newtonian (hereafter PN) equations of motion written in the center of mass frame including all 
terms up to 2.5PN order: 

^ = [(1 + A)n + 0v] + 0{l/c% (1) 

at 

where n = r/r, the coefficients A and B are complicated expressions of the binary's relative 
separation and velocity \27\ [28] . The Post-Newtonian approximation is a power series expansion 
in 1/c: the 0*'' order term corresponds to the dominant Newtonian acceleration; the even order 
IPN and 2PN terms are conservative and proportional to and respectively, and they 
are responsible for the precession of the pericenter; finally, the lowest order dissipative 2.5PN 
term which is proportional to 1/c^ causes the loss of orbital energy and of angular momentum 
due to radiation reaction. We treat the MBHs as point particles and thus we neglect any 
spin-orbit or spin-spin coupling. The runs were carried — either with the direct-summation 
NB0DY4 [29] or (^-GRAPE [30] codes — on the high-performance GRAPE-6A cluster at the 
Astronomisches Rechen-Institut (Heidelberg). The MBH pair motion, once it becomes a hard 
binary, is integrated with the 4*^^ order Hermite scheme, after performing a transformation to KS 
regularized variables [29[ 131] . The PN corrections are treated as perturbations to the dominant 
Newtonian acceleration in exactly the same manner as the perturbations from the field of stars. 
In Fig. 2, we show that both type of perturbations are, up to the very late stages of the 
relativistic inspiral, always much weaker than the dominant Newtonian acceleration. Therefore 
we are safe in adopting the linear approximation when adding both contributions to obtain the 
total perturbation. Although the IPN and 2PN terms are conservative and thus do not drive 
directly the relativistic inspiral, their strength is still, for most of the time, several orders of 
magnitude higher than that of the dissipative 2.5PN term. As a result, these conservative terms 
should be included in the calculation at all times in order to get an accurate evolution for the 
orbital elements during the whole inspiral and until they reach the LISA band [12\ [52] . 

We next describe a run that is typical of our sample. In this fiducial run, we set the black 
hole masses Mi and M2 of the two MBH particles to be 1% and 2.5% {q = 0.4) of the total 
mass in stars. We employ a total of = 128K particles to represent the total number of equal 
mass stars (half on each nucleus). Both nuclei have central densities with logarithmic slope 
7 = 1. Since LISA will be sensitive to the inspiral of MBH binaries of total mass M12 < IO^Mq, 
we have chosen to scale our model in such a way that the lighter MBH has Mi = W^Mq in 
physical units. As a result, one spatial A'-body unit corresponds to 5 pc and one iV-body time 
unit corresponds to ~ 1.67 x 10"^ yr. The resulting speed of light in our model units is c = 1022. 
We place initially the centers of the two nuclei at a separation of 20 model units, i.e. 100 pc. 
The half-mass radius of each nucleus is, for the fiducial run with 7 = 1, ri/2 ~ 2.41 in model 
units, or some 12 pc; this means that the nuclei are well separated in their initial configuration. 
The total mass in stars, initially enclosed within the inner 12 pc of each nucleus, is 5 x IO^Mq. 

After roughly 80 A-body time units, the MBHs have reached the center of the remnant and 
form a bound pair. At a first stage, they eject a large amount of mass (Mej ~ M12) in stars and 
its semimajor axis has shrunk by a factor of ~ 50 by time t = 100. Subsequently, the hardening 
rate slows down considerably but does never stall. 

The binary forms with an eccentricity close to 0.9, which is typical of our simulation sample, 
and grows slowly over time before it eventually starts to circularize due to radiation reaction. 
In Fig. 3a, we can see that the MBHs coalesce after roughly 450 A-body units, or some 7.5 
Myr. We should provide a word of caution, though, since the inspiral time scale is very sensitive 
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Figure 3. (a) Evolution of the inverse semimajor axis; (b) evolution of the eccentricity; (c) the 
evolution of the intermediate axis ratio h/a\ (d) the evolution of the minor axis ratio c/a. The 
different colors represent five different radii: r = 1 (red), r = 3 (green), r = 5 (dark blue), r = 7 
(magenta), and r = 9 (light blue). 

to the eccentricity evolution as a function of the binary's semimajor axis |26j . The study of 
the dependence of the evolution of the eccentricity as a function of the initial conditions, and 
especially as a function of the total number of stars (which, for a fixed total mass, determine 
the mass ratio between the MBHs and the stars) is the subject of ongoing work |22j . 

In panels 3c and 3d, we can observe the evolution of the axis ratio for the mass distribution 
of the merger remnant at several different radii. The radius of influence of the binary — 
radius enclosing 2Afi2 of the stellar mass — is in A'"-body (physical) units 0.2 (5 pc). The overall 
shape of the cluster is markedly triaxial throughout the simulation; at the inner radius r = 1 
represented in the figure, the triaxiality parameter T = (a^ — 6^)/(a^ — c^) decreases faster 
but is still ~ 0.1 — 0.25 by the time the binary decouples from the stellar environment at time 
t - 410 -420 (6.9 Myr). 

As the semimajor axis of the binary shrinks through the ejection of stars, the central density 
decreases progressively. If the triaxial nature of the potential drives the collisionless refilling of 
the loss cone, it is expected that a fair amount of stars supplied to the binary originate from 
a region well outside its influence radius rh |10] . On the contrary, in the case of an empty loss 
cone in a spherical symmetric nucleus, the scouring effects on the central density cusp would be 
limited roughly to the region interior to rh [33J. In our model, the influence radius is initially 
rh ~ 0.18 and evolves only very slowly to rh ~ 0.22 at the end. From the inspection of Fig. 
4, we can immediately see that the density cusp is being depleted out to radii of ~ fewxr/^, 
and this is a generic result in our simulation sample [22]. Furthermore, we measure the amount 
of mass ejected by the binary, by t = 320, from the region interior to Svh to be « 2M12; in a 
spherical nucleus, we expect a much smaller value M^j ~ 0.5Mi2 [33]. Note that this high rate 
of mass ejection may help to explain the existence of large mass deficits in some bright elliptical 






Figure 4. The scouring effect 
on tlie central density cusp of 
the merger remnant as the MBH 
binary semimajor shrinks through 
the ejection of stars. The mass loss 
extends to regions well outside of 
r/j, as expected if the loss cone is 
kept full by a centrophilic family of 
orbits associated with the triaxial 
stellar potential. The arrow signals 
the binary's radius of influence. 



galaxies. 

We can also observe that towards the end of the run, when the MBH pair has almost 
completely decoupled from the galactic nucleus (t ~ 410) and is evolving mainly in isolation 
under the effect of radiation reaction, the density cusp does not decay anymore. In fact, the 
cusp is starting to grow again, though at a much slower rate than that with which it decayed 
before, as expected in the case of a single MBH. 

Can we extrapolate these results to real galaxies? The answer to this question depends 
on the extent to which the evolution of the orbital elements (hardening rate and eccentricity 
evolution) converges for the range of particle particle number used in our simulations. Our 
preliminary results indicate that the hardening rate does indeed converge once we reach 
N ~ 0.125 — 0.25 X 10^. This is presumably due to the fact that the repopulation of the 
loss cone orbits happens in a regime that is essentially collisionless as a result of the nonlinear 
dynamics of a centrophilic family of orbits supported by the triaxial remnant nucleus. 



1e-12 
1e-14 


n=1 

n=2 

n=3 

'n=4 


Mbh=3.5x10^Mo ■ 
z=4 


1e-16 






„ 1e-18 






1e-20 






1e-22 






1e-24 






1e-26 







1e-07 1e-06 1e-05 1e-04 0.001 0.01 0.1 1 
f(Hz) 

Figure 5. Locus of the MBH binary, 
for the first six harmonics, in the LISA 
sensitivity band during their final inspiral and 
coalescence. The dimensionless strain he is 
plotted against the observed frequency. The 
total mass of the binary is Mi2 = 3.5xlO^M0. 
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Figure 6. The same as Fig. 5, but for 
Mi2 = 3.5 X IO^Mq. 



4. Gravitational Wave Signal 

The inspiral and merger of MBHs are expected to be one of the brightest sources of GWs 
to be detected by LISA [34j- In most of the literature, the strain amplitude of the GWs is 
estimated under the assumption of circular orbits. This has been motivated by the idea that 
MBH binaries should have been completely circularized by the time they enter the LISA band 
at / ~ 10"^ Hz [35]. In our runs, the binaries form typically with very large eccentricities and 
usually reach small separations (e.g. 100 Rschw) with non-negligible eccentricty. Depending on 
the redshift of the source, this implies that they may enter the LISA band with significant power 
distributed among the higher harmonics fn = ra/oj,b/(l + z), with n >2 (note that n = 2 for 
a circular orbit) [36J. In Figs. 5 and 6, we plot the characteristic strain amplitudes he for the 
our fiducial run when the model is scaled in such a way that the total mass of the BH binary is 
Mi2 = 3.5 X 10^(3.5 X 1O®)M0 and if the redshift of the source were 2 = 4. We conclude that 
several higher harmonics stand above the LISA sensitivity curve. We suggest that MBH binary 
orbits with non- vanishing eccentricity may be seriously considered for the LISA data analysis. 



5. Conclusions 

We present preliminary results from a set of TV-body models of merging galactic nuclei. The 
generic outcome of these mergers is the formation of a triaxial remnant. The MBHs present 
in these models often form highly eccentric binaries and, driven by the large pool of stars in 
centrophilic orbits associated to the triaxial potential, do not stall at a separation a ~ o/j. Those 
binaries that get formed with the highest eccentricity, e > 0.9, can decay from initial separations 
of ~ 100 pc down to gravitational-wave-induced coalescence within ~ 10^ yr. These binaries 
may reach the LISA band with non-negligible eccentricity and thus it may be important to 
consider finite-eccentricity effects in the data analysis. The amount of stellar mass ejected by 
the binary throughout the inspiral is higher than estimates made for spherical nuclei — this 
may help explaining the mass deficits observed in some of the brightest elliptical galaxies. 
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